Chaos in the Solar System 



1 



Chaos in the Solar System 

Myron Lecar, Fred A Franklin, Matthew J Holman 

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, Massachusetts 
02138; e-maih mlecar@cfa.harvard.edu 

Norman W Murray 

Canadian Institute for Theoretical Astrophysics, McLennan Physical Labs, University of 
Toronto, 60 St. George Street, Toronto, Ontario M5S 1A7, Canada; e-mail: 
murray@cita.utoronto.ca 

KEYWORDS: celestial mechanics, planetary dynamics, Kirkwood gaps 
ABSTRACT: 

The physical basis of chaos in the solar system is now better understood: in all cases investigated 
so far, chaotic orbits result from overlapping resonances. Perhaps the clearest examples are found 
in the asteroid belt. Overlapping resonances account for its Kirkwood gaps and were used to 
predict and find evidence for very narrow gaps in the outer belt. Further afield, about one 
new "short-period" comet is discovered each year. They are believed to come from the "Kuiper 
Belt" (at 40 AU or more) via chaotic orbits produced by mean-motion and secular resonances 
with Neptune. Finally, the planetary system itself is not immune from chaos. In the inner solar 
system, overlapping secular resonances have been identified as the possible source of chaos. For 
example. Mercury, in 10^^ years, may suffer a close encounter with Venus or plunge into the Sun. 
In the outer solar system, three-body resonances have been identified as a source of chaos, but 
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on an even longer time scale of 10^ times the age of the solar system. On the human time scale, 
the planets do follow their orbits in a stately procession, and we can predict their trajectories for 
hundreds of thousands of years. That is because the mavericks, with shorter instability times, 
have long since been ejected. The solar system is not stable; it is just old! 
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1 INTRODUCTION 



All known cases of chaos in the solar system are caused by overlapping resonances. 
Brian Marsden and I (ML), when we were students at Yale, had frequent discus- 
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sions with Dirk Brouwer about the significance of two resonances that overlapped. 

The "small divisors" that occur in perturbation theory at single resonances can 
be removed by a change of variables introduced by Poincaxe, which results in a 
resonance Hamiltonian similar to that of a pendulum. No such device was ever 
discovered for two overlapping resonances, and Brouwer sensed that there was 
something special about that case. Because a single trajectory can be numeri- 
cally integrated through an overlapping resonance without apparent catasrophe, I 
argued forcefully (and incorrectly) that the difficulty with overlapping resonances 
might just be a defect in the perturbation theory. However, had we numerically 
integrated a clone, initially differing infinitesimally from the original trajectory, 
the two would have separated exponentially. In fact, that is the definition of a 
chaotic orbit: exponential dependence on initial conditions. 

In many cases we can estimate the Lyapunov Time (the e-folding time in the 
above example) and even the Crossing Time (the time for a small body to develop 
enough eccentricity to cross the orbit of the perturber). Both times depend on 
the Stochasticity Parameter, which measures the extent of the resonance overlap. 

There is evidence that, today, small bodies in the solar system (e.g. comets 
and asteroids) behave chaotically. Meteorites are thought to be fragments of 
asteroid collisions. The asteroid Vesta has a reflection spectrum that resembles 
that of many meteorites. Every so often a meteorite hits the Earth, so we have 
evidence, within the last hundred years, of chaotic behavior. We have also found 
meteorites that originated on Mars. We believe they came from Mars because 
trapped bubbles of gas coincide with samples of the Martian atmosphere. 

Kirkwood (1867) noticed that the asteroid belt has "gaps" at resonances, i.e. at 
distances where the asteroidal periods are a rational fraction of Jupiter's period. 
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There have been many attempts to explain these gaps on the basis of a single 

resonance, but these attempts never produced gaps as devoid of bodies as the 

observed ones. We now understand that bodies are removed from regions where 

the overlap of two or more resonances induces chaos and large excursions in the 

eccentricity. Wisdom (1982, 1983, 1985) first illustrated how chaos at the 3:1 

resonance with Jupiter could result in sufficiently large eccentricity to allow an 

encounter with Mars or the Earth. Subsequent work showed that even collisions 

with Sun are a likely outcome (Ferraz-Mello &i Klafke 1991, Farinella et al 1994). 

Chaotic dynamics in the asteroid belt will be reviewed in Section 3. 

One to two "short period" comets (short means periods less than 200 years) 
are discovered per year. Short period comets are confined to the ecliptic and are 
believed to come from the Kuiper Belt, which is located about 40 AU from the 
Sun, in neighborhood of Pluto. Their stability is also discussed in this review. 

In contrast, long period comets are thought to come from the Oort Cloud at 
20,000 AU. They are perturbed into the inner solar system by stellar encounters or 
by the tidal field of the galaxy. These mechanisms differ from those determining 
the dynamics of short period comets and are not reviewed here. 

By now, the entire solar system exterior to Jupiter has been surveyed for sta- 
bility. Holman & Wisdom (1993) found that all the unpopulated regions of the 
solar system are unstable, on time scales that, in general, are much less than the 
age of the solar system. However, some comet orbits in the Kuiper Belt are just 
now becoming perturbed into the inner solar system. The chaotic dynamics of 
comets in the outer solar system is reviewed in Section 4. 

It is not too alarming that small bodies behave chaotically. Comets and as- 
teroids have individual masses less than 1/1000 that of the Earth and, in total. 
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make up much less than an Earth's mass. However, we depend on the regularity 

of the planetary orbits. Could they be chaotic? The answer for the known plan- 
ets is yes — but on a long time scale. The planets have presumably followed their 
present orbits for much of the lifetime of the solar system. But for how much 
longer? Are we in danger of losing a planet soon? Computers are just now able 
to integrate the planets for the life time of the solar system so we now have a 
preliminary exploration of this important problem. Here too chaos is induced at 
resonances, but evidence suggest that they are secular resonances that operate 
on long time scales. 

Earlier, the stability of the solar system was studied by looking for terms in the 
semimajor axis that grew with the time, (secular terms), or as the time multiplied 
by a periodic function of the time, (mixed secular terms). Now we know that 
instability comes from a chaotic growth of the eccentricity. The stability of the 
planets is reviewed in Section 4. 

For the reader who is intrigued, but new to this subject, we suggest a popular 
book called Newton's Clock: Chaos in the Solar System by mathematician Ivars 
Peterson (1993). We refer all readers to the earlier review by Duncan &; Quinn 
entitled "The Long-Term Dynamical Evolution of the Solar System," which ap- 
peared in the 1993 edition of this Annual Review. In addition, the proceedings of 
the 1996 workshop on "Chaos in Gravitational N-Body Systems" have appeared 
as a book and in Volume 64 of Celestial Mechanics and Dynamical Astronomy. 

Henon (1983) gives a lucid introduction to chaotic orbits and the "surface of 
section" technique. Ott (1993) has written the standard text on chaos which cov- 
ers a variety of physical problems. Sagdeev et al (1988) review other interesting 
applications of chaos including turbulence. 
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2 CHAOS AND CELESTIAL MECHANICS 

The rigorous condition for a mechanical system to be stable for all time is that 
there exists an "integral" (a conserved quantity) for each degree of freedom. The 
Sun and one planet is just such a system (e.g. the Kepler Problem). Fortunately, 
because the Sun is 1000 times more massive than Jupiter and the rest of the 
planets add up to less than a Jupiter mass, treating the planets as independent 
two-body problems is an excellent starting approximation. Using their "Keple- 
rian" orbits, one can calculate a first approximation to their forces on each other. 
Successive iterations of that procedure, carried out with great sophistication, are 
the techniques called celestial mechanics. The classic text on celestial mechan- 
ics was written by Brouwer & Clemence (1961). See Solar System Dynamics by 
Murray Sz Dermott (1999) for a recent treatment. 

Much of the history of celestial mechanics has involved the search for inte- 
grals of motion. The search was doomed to fail; eventually Poincare proved 
that there was no analytic integral for the problem of the Sun and two planets. 
However, there do exist nonanalytic integrals. Their discovery culminated in the 
Kolmogorov-Arnold-Moser (KAM) theorem, a fundamental result in the math- 
ematics of chaos. The theorem guarantees the existence of "invariant curves" 
(i.e. other integrals) as long as the perturbations are not too large and the cou- 
pling is not too near any resonance. Understanding the exact meaning of the 
word near was crucial, because resonances (like the rational numbers) are dense. 
This theorem is discussed in The Transition to Chaos by Reichl (1992) and in 
Regular and Chaotic Dynamics by Lichtenberg & Lieberman (1992) (our recom- 
mended text). For the more mathematically minded, there is also a set of lectures 
by Moser called Stable and Random Motions in Dynamical Systems (1973) and 
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Mathematical Methods of Classical Mechanics by Arnold (1978). 

Although the KAM theorem is of fundamental importance for the mathemat- 
ical structure of chaos, the strict conditions of the theorem are satisfied in the 
solar system. In what follows, we will be concerned with orbits that are not 
covered by the KAM theorem. 

All of the analytic work described in this review relies on perturbation theory, 
exploiting the fact that the planetary masses are small compared with the mass 
of the Sun. A brief outline of the approach follows. 

The equation of motion of a planet orbiting a star accompanied by a second 
planet is 

^ + g(Mo + Mi)^ = Vii?i,2, (1) 

where 

" 1 ri • r2 



Rl,2 = QM2 



(2) 



/i,2 n 

Here ri is the distance between the Sun and planet 1, and ri,2 = -\/(ri — r2)^ 
is the distance between the planets. The quantity R, which is the negative of 
the planetary potential, is known as the disturbing function; it describes the 
disturbances of the planet's elliptical orbit produced by the other planet. The 
second term in the square brackets arises from the noninertial nature of the 
coordinate system employed and is known as the "indirect" term. It occurs 
because the traditional coordinate system takes the position of the sun as the 
origin. It is not generally relevant to the chaotic behavior of bodies in the solar 
system, so we ignore it henceforth. 

The next step is to expand the disturbing function using the expressions for 
ri and 9i found by solving the Kepler problem. This rather daunting task has 
been performed by a number of authors (Peirce 1849, Le Verrier 1855, Murray & 
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Harper 1993). The general form is 

X cos [71 Al + j2A2 + jsZUl + ^4^2 + j5^1 + j6^2] ■ (3) 

In this expression the ji are positive and negative integers. The angles in the 
argument of the cosine are the mean, apsidal, and nodal longitudes and are 
measured from the x-axis. If we rotate the coordinate system, the disturbing 
function, which is proportional to a physical quantity (the force), cannot change. 
This implies that J2iji = 0- We have kept only the lowest order terms in the 
sum; for a given terms proportional to e'-'''"'"^ or larger powers will also appear. 
Each cosine represents a resonance; the effects of these resonances constitute the 
subject of this review. 

We can use this expression to find the effect of one planet on another. Hamil- 
ton's formulation of mechanics offers the easiest way to proceed. Using action- 
angle variables for the two body problem, the quantities A, —w, and — J7 are 
appropriate angles. The corresponding actions are simple functions of a, e, and i: 
L = ^/gM^, G = ^/g^^L[l-^/{l-e^)] « {l/2)e'^L, and H = ^/gMa{l - e^){l- 
cos i) « (l/2)i^L, where we have distinguished the gravitational constant, Q, from 
the momentum, G. In these variables the Kepler Hamiltonian isH = {QM)/2L'^ . 
Because we are using action-angle variables, none of the angles appear. This 
tells us immediately that the motion takes place on three-dimensional surfaces 
in phase space, defined by L, G and H held constant. Topologically, this is a 
three-torus. The surprising thing is that G and H do not appear in the Hamil- 
tonian, so that w and J7 are also constant. In a generic three degree of freedom 
Hamiltonian one would expect all three actions to appear explicitly, leading to 
three non-constant angles. The Kepler problem is degenerate, as is reflected in 
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the fact that H and G do not appear in 7i. In this case the motion takes place 

on a one-torus, or circle, in phase space. 

To find the variation in a, for example, note that 

dLi 1 Li dai 



dt 2 ai dt ^ 



(4) 



so 



doi _ (K) 
dt ~ LidXi' ^ ^ 

The presence of planet 2 forces periodic variations in the semimajor axis of planet 
1. For example, suppose we pick a term of the form 

</'2 -5,3,o,o,o(ai, 02)6? cos[2Ai - 5A2 + 3tz7i]. (6) 

The equation for oi becomes 

^ = _4(a^ni) — ^(/)2 -5,3,0,0,0(01, a2)ef sin[2Ai - 5A2 + 3zui] + 0(e)^. (7) 
at 02 IVl 

Similar expressions can be derived for all the other orbital elements. 

With these expressions, or similar ones provided by Lagrange, it appeared to 
be a simple matter to integrate the equations of motion. However, early efforts 
to do so revealed difficulties. The problem can be seen by integrating equation 
d?!). To first order in M2 we find 

ai{t)-ai{0) = ai — ^- , ^. '/'2,-5,3,o,o,o(qi, Q2)ef cos[2Ai-5A2+3CTi]+0(e)^. 

(8) 

The denominator 2ni — 5n2 + 3wi, and similar denominators that arise when all 
the terms in the disturbing function are considered, is the source of the difficulty. 
Poincare (1993) pointed out that for integers ji and j2 (2 and 5 in our example) 
the denominator becomes arbitrarily small. He went on to show that, in spite 
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of the fact that terms with large ji and j2 tend to carry large powers of the 

eccentricities, the sums used to define the hoped-for analytic solutions diverged. 

The locations in phase space (essentially along the a axis) where the denomi- 
nator vanishes are known as resonances. Poincare noted that in the immediate 
vicinity of a resonance the motion was very complicated. Later work, particu- 
larly that of Chirikov (1979), showed that the complicated motion, dubbed chaos, 
occupied a large fraction of the phase space near the resonance when two neigh- 
boring resonances overlapped. This result is essentially an anti-KAM theorem, 
in that it specifies where KAM tori are absent. It provides the underpinings of 
much of the work reviewed in this article. 

2. 1 The Origin of Chaos: Overlapping "Mean-Motion" Resonances 

The chaos generated at overlapping resonances was first studied, in the astronom- 
ical context, by Wisdom (1980), who calculated the width and extent of overlap 
of adjacent first-order resonances in the circular restricted three-body problem 
and found that they overlap to a distance given by (5a/a ~ l.S^u^/^, where is the 
mass ratio of the planet to the central star. The resonance overlap criterion for 
chaos was developed by Chirikov (1979). The width of the first-order resonances, 
for zero eccentricity, in the restricted problem was also derived by Franklin et 
al (1984), from which we find that the 2:1 mean-motion resonance extends, in 
semimajor axis (with ajupUer = 1) from 0.621 to 0.639, the 3:2 from 0.749 to 
0.778, the 4:3 from 0.808 to 0.843, the 5:4 from 0.841 to 0.883, and the 6:5 from 
0.863 to 0.908. The 4:3 overlaps with the 5:4 and all adjacent resonances closer 
to Jupiter overlap. The 4:3 resonance is at 0.825; the Wisdom formula predicts 
onset of chaos for a > 0.82. This region where mean-motion resonances overlap 
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is where the relation between the Lyapunov time, Tg, and the time for a close 

encounter with the perturber, Tc, applies. This relation, found empirically by 

Lecar et al (1992), predicts that Tc, is proportional to 7]^'^^. 

Murray &; Holman (1997) have explained this relation in terms of the Stochasity 
Parameter, K = (vrAA/M)^, where AA is the width of the resonance and SA 
is the separation between the resonances. If AK = K — Kc, where Kc is the 
critical value (~ 1, corresponding to the start of overlap), then they showed that 
Ti ~ AK-^-^^ and ~ AiC-^-^s, so ~ T/^ ^ This holds in the region of 
overlapping first-order resonances (the /x^/^ region) . Orbits in this region show a 
range of more than three orders of magnitude in Tg. However, outside that region, 
they showed that the relation was in error by a factor of 10 at the 5:3 resonance 
(a second-order resonance), and by a factor of 100 at the 7:4 resonance (a third- 
order resonance). They also integrated 10 "clones" of Helga, an asteroid at the 
12:7 resonance (a fifth-order resonance). Five of the clones had Lyapunov times 
ranging from 6,000 to 13,000 years. They encountered Jupiter in 1-4 Gyr (Gyr = 
10^ years). In this case the relation predicted Tc ~ 6 Myr (Myr = 10^ years) — too 
low by a factor of 1,000. Murray &: Holman (1997) showed that the mechanism 
for chaos in these higher order resonances was an overlap of the "subresonances," 
and that diffusion between overlapping subresonances in the same mean-motion 
resonance is slower than diffusion between overlapping mean-motion resonances. 

It is worth noting that higher order resonances become very narrow for zero 
eccentricity. For example, from Franklin et al (1984), the width of first-order 
resonances is proportional to fi^^^ for zero eccentricity. The corresponding widths 
of second-order resonances is proportional to /x and the width of third-order 
resonances is proportional to fi^. Thus, second and higher order resonances are 
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too narrow to overlap each other at low eccentricity. However, higher order 

resonances can occur in the wings of first-order resonances. 



3 DYNAMICS IN THE ASTEROID BELT 

Understanding the distribution of the asteroids as a function of their semima- 
jor axes, and particularly where their mean-motions are commensurate with 
Jupiter's, provided dynamicists an intriguing puzzle for over 130 years — all the 
more so because these (Kirkwood) gaps occur at most mean-motion ratios, UA/njup = 
p/q, but a concentration of bodies at two others. Progress on this classic problem 
has been striking over the past 20 years. In a broad sense, it has been solved: We 
can identify the sources of orbital instability (or their absence) and the nature of 
their consequences and also have a good idea of some of the time scales involved. 
Numerical and analytic studies both have contributed extensively. Although sev- 
eral related dynamical processes have been — and still are — working to produce 
gaps in the asteroid distribution, the most significant ones can all be linked to 
the solar system's present environment. Carving gaps may in some cases require 
upwards of a billion years, but it can probably be done without requiring cos- 
mogonic explanations; i.e. calling on processes that occured in the primordial or 
developing solar system. 

The paper that ignited the modern era of work on the Kirkwood problem was 
Jack Wisdom's (1982) first contribution to the study of the 3:1 mean-motion reso- 
nance at a = 2.50 AU. His startling results showed that an orbit at this resonance 
could remain quiescent, with a low eccentricity, e < 0.1, for more than 100,000 
years but also show occasional surges lasting about 10,000 years that would lift 
e to a maximum value of about 0.35. Such a value is just sufficient to allow a 
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crossing of Mars' orbit, resulting in an eventual collision or a close encounter. 

Orbital computations as long as a million years were a rarity 20 years ago, and 
Wisdom's novel approach was to develop a mapping of the planar elliptic three- 
body problem that relied on two efficient techniques. The first approximated the 
short-period terms (i.e. ones characteristically arising during an orbital period) 
by a series of delta functions. The second averaged the Hamiltonian, expanded to 
second order in eccentricity, over the longer but still relatively short-term angular 
variable that librates at the 3:1 resonance. 

In two following papers. Wisdom (1983, 1985) first used direct numerical in- 
tegrations to verify the presence of the emax peaks of 0.35 and, by including a 
third dimension, then showed that Cmax could rise to e = 0.6, a value that also 
included Earth-crossing trajectories. At the same time he calculated the extent 
of the chaotic zone, showing that it closely matched the observed 3:1 gap width. 
The excitement generated by these results echoed widely: A straightforward dy- 
namical process that (a) could open a gap at one resonance, (b) in principle might 
be generalized to account for other Kirkwood gaps, and as an added bonus, (c) 
could deliver asteroidal fragments into the inner solar system as meteorites, had 
at long last been identified. Later work by Ferraz-Mello & Klafke (1991) and 
Farinella et al (1994) showed that the chaotic zone at low e is linked, even in the 
planar elliptic three-body problem, to one with e > 0.6 so that e — > 1 can occur. 

The panels of Figure || provide a quick insight into the chaotic behavior that 
Wisdom discovered at 3:1. Figure |l|a plots the motion of both Jupiter's apsidal 
line, zuj, and that of a low eccentricity body [cq = 0.05] in the 3:1 resonance, 
and Figure shows that the eccentricity surges occur when the two ro's are 
approximately equal — in fact, vja > t^j corresponds to the rise in e. The equality 
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of two apsidal or nodal rates is referred to in solar system studies as an example of 

secular resonance. Their importance at certain locations in the asteroid belt has 

long been recognized (cf Brouwer & Clemence 1961), but only in the past decade 

has their role within mean-motion resonances been appreciated. Perturbations 

arising at mean-motion resonances will markedly effect the elements of bodies 

lying in them — yielding a very broad range of apsidal and nodal rates that are 

functions of e and i. The apsidal and nodal motion of the Jupiter-Saturn system 

is defined to high precision by two apsidal and one nodal terms that have been 

labeled 1/5, uq, and i^ig by Williams (1969). Figure |^ shows that the apsidal 

motion of a body with ao = 0.481 and cq = 0.05, which librates in the 3:1 

resonance, will intermittently also resonate with the frequency of 1^5 and possibly 

as well. 

The challenge to map the locations and limits of the secular resonances that lie 
within the confines of many significant mean-motion resonances has been met on 
theoretical grounds in papers by Moons & Morbidelli (1995 and its references). 
Although this work concentrates on the planar case in which only z^s and uq are 
present, it contains the important result that, in the higher order resonances, 
3:1, 5:2, and 7:3, the z/5 and secular resonances exist over a very wide and 
overlapping range of a and e — to such an extent that a condition of widespread 
chaos is present inside these three mean- motion resonances. Figure^ is a sample 
of their work for the 5:2 resonance. 

The picture developed by Moons & Morbidelli (1995) provided one reason that 
led Gladman et al (1997) to study the fates of a large number of bodies placed 
in various mean- motion resonances. In a real sense their paper represents an 
elaboration and even a culmination of Wisdom's original suggestion that marked 
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eccentricity increases are responsible for the Kirkwood gaps. Some eccentricity 

increases beyond the value of 0.35 found by Wisdom had already been noted, 

but Gladman et al provided accurate statistics by integrating more than 1000 

bodies sprinkled throughout 3:1, about 450 in 5:2, and 150 each in 7:3, 8:3, and 

9:4. Their study quantitatively describes the dynamical transfer process noted 

earlier by Wisdom (1983) — namely that gravitational encounters even with the 

terrestrial planets can provide sufficient energy changes to move bodies from 

regular to chaotic zones and even from one resonance to another. We can now 

legitimately claim that the development of a gap [cf Figure ^ at 3:1 is inevitable, 

though some details are complex and different time scales are followed. The next 

two paragraphs provide an outline of the process. 

First, for bodies once in 3:1 with e < 0.25, the effect of imbedded secular 
resonances, principally uq, will drive e's of any and all bodies toward unity, leading 
most likely to solar impacts in times of a few million years. In the survey, this 
was the fate of about 70% of the initial population. A quarter to a third of them 
directly impacted the Sun, whereas the majority were gravitationally scattered 
by the Earth or Venus before doing so. Most of the remaining 30% moved in 
unstable orbits exterior to Saturn. 

Second, the same fate awaits bodies of any eccentricity whose critical arguments 
are either circulating or that have librations greater than about 50°. However, 
bodies with e > 0.3 and that show small librations have different outcomes. In 
preparing for this review, we identified orbits of bodies with 0.3 < e < 0.6, whose 
moderate librations, all < 40° qualified them as stable librating members of the 
3:1 resonance. These orbits are regular; i.e. there is no sign of any exponential 
growth in their angular orbital elements during integrations lasting as long as 10^ 
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yrs. Figure]^ — c is an example of one of them. Three (of three) remained in reg- 
ular orbits until the integrations were terminated after 2 billion years. However, 
the eccentricities of all three of these "stable" orbits regularly climbed at least to 
0.5, hence risking the gravitational encounters (or, less probably, actual collisions) 
with Mars mentioned above. (Mars itself was not included in these integrations.) 
Gladman et al found that 5% of all bodies initially in 3:1 were "extracted by 
Mars," meaning that their orbits were first perturbed by Mars by a sufficient 
amount that the final result after subsequent encounters was most likely a solar 
impact or an orbit beyond Saturn. This seems the certain fate of the otherwise 
stable high eccentricity bodies — "otherwise" meaning the case with only Jupiter 
and Saturn present. Their result argues that this phase of the depopulation of 
3:1 will require longer times, 10-100 Myr (with a tiny handful remaining after 
100 Myr, but the eventual outcome is the same as that of their lower eccentricity 
neighbors. We can conclude that 3:1 is a resonance that has been emptied of any 
asteroids initially present by a natural, multi-stage dynamical process in which 
all planets, Venus through Saturn, have contributed. 

The 5:2 resonance, although one order weaker than 3:1, behaves similarly but 
with some interesting differences. As was the case at 3:1, values of the eccentricity 
set three regimes: (a) Orbits of low e are severely chaotic owing to the influence 
of fg. This remark applies to all bodies with a ~ 05:2 = 2.78 AU and all e < 0.2. 
Figure ^a-d presents an example, {b) for 0.2 < e < 0.40, orbits are regular, 
provided that their libration amplitudes are less than 11° and (c) for all larger 
e's, orbits at least occasionally lie in secular resonance and are very chaotic with 
e's reaching values > 0.7. These remarks emphasize the resemblence to 3:1, but 
there are two novel features for bodies in the second category. 
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First, their e ranges generally do not include a close approach to Mars and 

second, as figures show, in just one of many cases, the presence of secu- 
lar resonance does not always correspond to extreme chaos. Taken together, 
figures |5|a-c? are examples of the differing response of two bodies at the 5:2 mean- 
motion resonance to secular resonances. When the amplitude of ro^i — luj is large 
or when there is coincidence of the periods of the term with the uja oscilla- 
tions, then chaos is severe and is measured by Lyapunov times of ~ 1000 Jovian 
periods. By contrast, when the oscillations of tua no longer match the frequency 
of and/or their amplitude is small the orbits remain regular. A tentative con- 
clusion is that z^5 is far weaker than fg. Thanks to Moons & Morbidelli (1995) 
we have an accurate knowledge of the locations of the secular resonances as well 
as the limits on the orbital elements over which they operate. What we lack is 
an evaluation of their strength — say, a measure of the chaos they can produce 
within a mean-motion resonance. Murray & Holman (1997) have developed the 
means to calculate Lyapunov times when mean- motion resonances overlap. A re- 
lated formalism applying to the case when secular and mean-motion resonances 
overlap would be very valuable. 

In their study Gladman et al found that the loss of bodies from 5:2 closely 
resembled the depletion at 3:1; i.e. all objects were removed within 100 Myr. 
Our suggestion is otherwise: We expect a number of librating bodies with e < 
0.40 to remain, unless a slow diffusion into the chaotic region has removed them 
all. However, there may not be a real conflict. Gladman et al were especially 
interested in the role of the 5:2 in delivering members of three nearby asteroid 
families into the inner solar system and hence chose initial velocity distributions 
accordingly. Inclusion of librators in such a distribution is very unlikely. By 
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contrast, their survey at 3:1, though it also followed the evolution of members of 

three families, introduced 1000 randomly selected orbits as well. The files of the 

Minor Center list some 30 possible candidates in the above eccentricity range at 

5:2. Whether they do librate with small amplitudes and avoid secular effects is 

at present unknown. 

Moons & Morbidelli (1995) have also mapped the one (weaker) secular reso- 
nance, f5, that lies within the 4:1 mean- motion resonance, but 4:1 itself is unim- 
portant because its semimajor happens nearly to coincide with the locations of 
the far stronger uq and i^ig secular resonances that, independently of any mean- 
motion resonance, are now known to be the agents defining the inner edge of the 
asteroid belt at 2.1 AU. 

Both Moons & Morbidelli (1995) and Gladman et al (1997) also examined the 
7:3 mean-motion resonance. The former paper again shows that both z^s and vq 
are centrally situated within 7:3 and cover large portions of it, both as e — > and 
at e as high as 0.65 as well. In preparing for this review, we randomly introduced 
120 bodies into the resonance and found severe chaos {Ti < 2000Pj) everywhere. 
Gladman et al (1997) found that depletion at 7:3 proceeds at a slower rate than 
at 5:2, with nearly one half of its initial bodies surviving to their integration limit 
of 40 Myr. Our results are roughly compatible with theirs but suggest a depletion 
rate that is faster by a factor of at least two. In any event, as Figure ^ attests 
and both results predict, 7:3 should be a genuine, fully developed gap, not just a 
region of reduced population. 

The fifth-order resonances, 8:3 (a = 2.71 AU) and 9:4 (a = 3.03 AU) were not 
examined in the otherwise comprehensive studies of Moons & Morbidelli so we 
have introduced about 50 bodies in each to have some idea where chaos is most 
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important. The results contain few surprises: at 8:3 orbits with Cp < 0.19 are 

severely chaotic, with T/'s < 2000Pj. In the range 0.19 < Cp < 0.35 they are quite 

regular, having logTi > 4.5 (Jovian periods). They become increasingly chaotic 

at higher e's, a fact that is of little importance because e = 0.4 leads to a Mars 

encounter. At 9:4 orbits are again very chaotic for Cp < 0.2 and much less so for 

0.2 < Cp < 0.32. Still, their average T;'s are shorter by a factor of about three 

than those at 8:3. Eccentricities > 0.32 are again very chaotic, and those with 

e > 0.45 will be Mars crossers. Gladman et al (1997) found incomplete depletion 

at both 8:3 and 9:4, with about one half of the original population surviving after 

40 Myr in the former and three quarters after 120 Myr at the latter. All of these 

estimates are qualitatively in accord with the distribution shown in Figure ^. [A 

likely interaction between 9:4 and 11:5 (3.08 AU) probably accounts for the extra 

width.] Longer surveys, with an eye to further quantifying diffusion rates into 

chaotic areas, are an important future project. 

3. 1 Chaos in the Outer Belt 

The 2:1 resonance at 3.28 AU divides the populous inner belt from the much 
less dense outer portion [cf Figure |3|. Can it be that the entire outer belt is 
systematically more chaotic than the inner belt? The answer is certainly yes, 
and it has recently become possible to evaluate both its extent and severity 
throughout all of the outer belt where Holman & Murray (1996) have, among 
other results, shown that 22 of 25 outer belt minor planets have Lyapunov times 
shorter than 6OOO-P7. However, it is also important to stress that there is more 
chaos in the inner belt than what we have thus far mentioned. Independent 
studies by Murray et al (1998) and Nesvorny & Morbidelli (1998) point out 



20 Lecar,Franklin,Holman &:Murray 

that many three-body mean- motion resonances (cf Aksnes 1988), involving the 

longitudes of Saturn as well as Jupiter and the asteroid, are a major source of 

chaos. Nesvorny and Morbidelli (1998) trace the chaos in the orbits of about 250 

of just the numbered minor planets to three-body resonances, finding Lyapunov 

times that lie in the range of 1,000 to 10,000 Pj. As the authors demonstrate, 

these resonances are vastly more dense in the outer belt, where their effect can 

only be more destabilizing. 

Holman & Murray (1996) and Murray & Holman (1997) have studied chaos in 

the outer belt by examining orbital behavior at a number of high-order resonances 

(e.g. 12:7) where they have calculated both Lyapunov and diffusion time scales, 

analytically as well as numerically, and some escape times for comparison. They 

associate chaos with the overlap between members of individual mean-motion 

resonances. Consider the example of 12:7. In the planar case it is composed of 

six (sub)resonances, with various multiples of the two apsidal longitudes, e.g. 

(Tfc = 12Aj - 7\a - kwA + {k- 5)wj, (9) 

for k = 0,1, ...5 and where the A's and tzj's refer to mean and apsidal longitudes. 
Two components, defined e.g. by k and k-|-l can overlap and Holman & Murray 
(1996), having obtained their locations and widths, showed that their overlap 
will generate chaos. This approach to analyzing chaos within a given mean- 
motion resonance is the mathematical equivalent of the one used by Moons & 
Morbidelli (1995), who consider the overlap between a resonance defined by a{k) 
in Equation |9| and the secular one given by 

CTsecular = WA-^J = [(y{k + 1) - CF{k)]. (10) 

The latter is the approach taken by Wisdom in his analysis of the 3:1 resonance 
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(Wisdom 1985). In this approach it is sometimes assumed that there is an ad- 
equate separation between the time scales associated with cr^k) and a{secular) 
so that the action associated with a{k) is adiabatically preserved, except when 
the orbit is near the separatrix. However, these time scales are not adequately 
separated for most resonances in the asteroid belt. 

Figure ^ plots some of Holman & Murray's results, indicating that (a) asteroids 
in the outer belt have already been ejected from resonances of order less than 4 
(with the exception of 11:7), thus demonstrating that less well-defined Kirkwood 
gaps will also exist there; (6) escape times from fifth-order resonances correspond 
roughly to the solar system's age' and (c) objects in sixth-order resonances are 
ejected in times as long as 10^^^^^ years. We encounter again the problem dis- 
cussed elsewhere in this review: Lyapunov times are measured in a few thousand 
Pj, but at least some of the carefully evaluated escape times exceed the age of 
the solar system. In view of these long times for ejection, it seems likely that 
other processes have contributed to the removal of some resonant and many non- 
resonant bodies from the extended outer belt. Holman & Murray (1996), Liou & 
Malhotra (1997), and Nesvorny & Ferraz-Mello (1997) have directed attention to 
planetary migration as a mechanism that will move various resonances into and 
through the outer belt, a process that can have further dynamical consequences, 
as the latter paper stresses, if there is also a change in the near 5:2 commensura- 
bility in the mean motions of Jupiter and Saturn. Lecar &: Franklin (1997) and 
Franklin & Lecar (2000) have quantitatively studied the sweeping of secular res- 
onances associated with the decay of the solar nebula through the asteroid belt. 
These papers show that this one mechanism can accomplish three desired ends: 
(a) remove an overwhelming fraction of an initial population from the outer belt. 
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(b) deplete the inner belt from a likely early value by a factor of about 1000 so as 

to match present observations, and (c) generate a range of eccentricities that are 

characteristic of the known minor planets. Their study also included gas drag on 

the asteroids. 

3.2 Behavior at First- Order Mean-Motion Resonances 

We turn at last to the behavior at two first-order resonances that Morbidelli &: 
Moons once called the "most mysterious ones," beginning with a discussion of 
the reason(s) for the pronounced gap at the 2:1 resonance, oq = 0.630 and the 
concentration of minor planets at the 3:2 resonance, oq = 0.763. Part but not 
all of this mystery has been dispelled as models have become more realistic. To 
be more precise, unless models include Jupiter and Saturn in eccentric precessing 
ellipses, they can not capture enough physics to account for the behavior at 
the first-order mean-motion resonances. The importance of secular resonance in 
developing chaos within such higher order mean- motion resonances as 3:1 and 
5:2 might lead to the expectation that they should also be of prime importance 
here as well. However, a detailed mapping (Morbidelli & Moons 1993) shows 
that the role of ly^ and z^g for generating chaos in the heart of both the 2:1 and 
3:2 resonances applies only for orbits having eccentricities e > 0.45 and 0.25, 
respectively. The reason is clear enough: At low e's, first order resonances drive 
pericentric longitudes, zua, much more rapidly than the motion of even the faster 
of the two principal terms (i.e. z^q) that measure vcj — by a factor of about 30 for 
orbits with e's averaging about 0.15. We must look elsewhere for other resonances 
whose overlap at orbits of low e with the 2:1 resonance itself is the source of the 
chaos shown in figure |^a. 
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The ones we are looking for involve commensurabilities between multiples of 

the (resonant) libration and apsidal frequencies. The former are defined by the 
regular oscillations of the appropriate critical angle a similar to Equation ^. 
Strictly periodic solutions have o" = 0; real bodies at mean-motion resonances 
show oscillations in a with well-defined periods and amplitudes that also corre- 
spond to regular variations in their orbital elements, especially semimajor axis 
and eccentricity. Periods of a increase with e, and at 2:1 and 3:2, typically lie in 
the range of 15 to 75 Pj, where Pj, the orbital period of Jupiter, is 11.86 years. 
Apsidal periods also depend on e, from close to 15Pj for very small values, then 
increasing rapidly. Therefore, a series of commensurabilities, starting near unity, 
which we label Papse/Pub in Figures |^ and |8| must occur. 

Giffen (1973) first called attention to these commensurabilities, though their 
link to chaos was not stressed. Henrard and colleagues (cf Henrard and Lemaitre 
1986, Lemaitre & Henrard 1990 first recognized their importance and began to 
investigate them analytically, and Franklin (1994, 1996), numerically. Following 
Henrard we shall refer to the small integer commensurabilities of Papse/Pm as 
secondary resonances. We can explore their behavior and importance by consid- 
ering the planar example shown in Figure |^a-c. The figure shows that they are 
indeed the source of chaos within the 2:1 resonance by presenting three cases in 
which the two perturbers, Jupiter and Saturn, move in planar precessing orbits 
with (case a) their present eccentricies (ej = 0.044 it 0.016; es = 0.047 ± 0.035), 
then with one half (case b) and finally one tenth (case c) of these values. This se- 
quence helps verify the claim that secondary resonances are responsible for chaos 
among (hypothetical) resonant asteroidal orbits of low eccentricity. [A slash, e.g. 
5/3, denotes a secondary resonance, while a colon, e.g. 3:2, denotes a mean- 
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motion resonance.] Figure shows that only the strongest secondary resonances 

are present and their widths are very narrow, but in figure when e(J) and e(S) 

have normal values, higher order secondary resonances [e.g. 5/3] have also been 

excited and/or first-order ones have broadened so that extensive overlap occurs. 

A quick interpretation of figure |7| might suggest that secondary resonances are 
the source of chaos in the entire range, < e < 0.28. If this is truly the case, then 
ones as high as 15/1 must still be contributing. Despite our present ignorance 
of their strengths, this in itself seems curious, though since Papse is <^ rapidly 
rising function of e(asteroid), whereas Pnh is almost constant, it is also true that 
higher order resonances crowd closer and closer together. The fact that exactly at 
some weaker secondary resonances very regular orbits exist is also curious, though 
helpful, as it speeds the process of locating them. A supplementary explanation 
for chaos near e = 0.2 has been mentioned by Morbidelli & Moons (1993), who 
found that the z/ig secular resonance is present in the limit as z ^ and should 
affect orbits with e = 0.21 + 0.04/ — 0.02. However, the calculations leading to 
figure 1^ are strictly limited to the planar approximation so that the origin of 
chaos for cq < 0.15 is clearly the province of secondary resonances, but to what 
e they extend is less clear. 

Recently, Moons et al (1998) broadened this topic by mapping the locations of 
chaotic zones that arise from all secular and families of secondary resonances at 
inclinations of 0, 10, 20, and 30°. [As Henrard (1990) noted, secondary resonances 
are not confined to the case defined by ci = wa^ but can include linear combina- 
tions of WA with the frequencies of planetary apsidal and nodal motion.] What 
is striking is the connected nature of the chaotic zones, extending through much 
of the area in the a, e plane for all four z's. Depending on the inclination, one to 
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three islands are present that contain quite regular orbits, having 0.2 < e < 0.45, 

many of which are not at risk because their perihelia avoid crossing Mars' orbit. 
Moons et al note that one of them is populated by five recently discovered minor 
planets, but the others are seemingly empty. 

With the source of chaos identified, the question now turns to how the deple- 
tion of a hypothetical early population of bodies at 2:1 might have proceeded. 
This topic rests on less firm ground. Two independent sets of long-term inte- 
grations argue that objects with Lyapunov times as short as several hundred 
Jovian periods (a few thousand years) cannot permanently exist at 2:1. Franklin 
(1996) found examples of escape for single representative bodies placed in the 
2/1, 3/1 and 5/1 secondary resonances, all after times close to 800 Myr. Objects 
in much higher secondary resonances, consequently with Lyapunov times three 
orders longer, remained after 4 Gyr. These integrations included only Jupiter and 
Saturn, which moved in planar precessing orbits. In five cases of escape orbital 
eccentricities rose to values that guaranteed a crossing of Mars' orbit only a few 
million years prior to escape, leading to the conclusion that drifting out of reso- 
nance, not an eventual encounter with Mars, is responsible for depopulating the 
region. At the same time Morbidelli (1996) reached a somewhat different conclu- 
sion, though one still compatible with the existence of a gap. His 3-dimensional 
integrations of 10 orbits with 0.055 < e < 0.155 indicated that all became Mars 
crossers in times between 10 and 100 Myr. (His integrations were terminated 
once this orbit crossing was noted.) The same fate happened to most orbits with 
higher e's, though four remained for a full 10^ years. Most objects remaining at 
that time showed strong evidence that their proper eccentricities were diffusing 
toward higher values. 
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The more rapid crossing times found by Morbidelli may be the consequence of 

adding the third dimension to the problem. In a parallel development Henrard 
et al (1995) showed analytically that a strong (nodal) secular resonance overlaps 
several secondary resonances. They suggested that it may be possible for a body 
initially in a low e secondary resonance to diffuse or random walk with an in- 
creasing libration amplitude so as to enter the extensive field of principally nodal 
secular resonances that lie at much larger e's and z's. Their short integration 
times of only 1 Myr failed to provide an example. At this point in our discussion 
the 1998 paper of Moons et al assumes special importance because it indicates 
that the chaotic areas arising from various sources are not isolated islands, but 
are connected. Bodies may therefore diffuse from low to high eccentricity and 
escape with only small changes in semimajor axis, much as was found in the 
numerical study just mentioned. How and when these three islands, located by 
Moons et al (1998) and containing very regular orbits, have lost bodies is, de- 
spite Morbidelli's long integrations, still an unsolved problem. We would argue 
for additional very long term studies to quantify the process of diffusion especially 
for orbits with small-to-moderate librations in the 0.2 < e < 0.45 range before 
placing all hope in an explanation linked to planetary migration (cf Nesvorny 
& Ferraz-Mello 1997). Despite our present state of uncertainty, there is now a 
general conviction that strictly dynamical effects linked to Jupiter and Saturn in 
their present orbits are still the best bet to exhaust a primordial population of 
bodies at the 2:1 resonance, but that such processes needed more than a billion 
years to yield the current minor planet distribution. 

Figures §a,6,c, the companion to Figures ^a, 6,c, present parallel results for 
the 3:2 mean-motion resonance, which is characterized not by a broad gap but 
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by a concentration of nearly 200 bodies (called the Hilda group) in reliable (i.e. 

observed for two or more oppositions) librating orbits. The difference in the 
observed appearance at the two resonances is clearly reflected in their dynamical 
behavior. Although the strongest secondary resonance at Papse/Pub = 2/1 is 
easily identified in Figure ^, it is less deep (less chaotic) and far narrower than is 
the case at 2:1. Even when the eccentricities of Jupiter and Saturn are doubled 
(cf.Figure ^ c) its width is far less than its counterpart at 2:1. A notable feature of 
the Hildas is the absence of bodies with proper eccentricities, Cp ^ 0.10. Figure 
shows that chaotic orbits dominate in the region eo — Cp < 0.05 and that a mix 
of chaotic and regular ones lie from eo = 0.05 to the 2/1 secondary resonance at 
eo = 0.074 with regular ones alone present at higher eo's, except exactly at other 
secondary resonances. 

Nesvorny & Ferraz-Mello (1997) have compared the long-term behavior through- 
out the 2:1 and 3:2 resonances. Their approach applies the frequency map tech- 
nique discussed by Laskar (1993) to determine the diffusion of (Fourier) frequency 
components in a resonant body's pericentric motion. The results show that all 
orbits in the heart of the 3:2 with 0.05 < e < 0.35 diffuse very slowly, with 
percentage frequency shifts of about 10% in a billion years. At 2:1 the broad, 
continuous region found at 3:2 is replaced by numerous smaller islands also hav- 
ing slow diffusion, but they are surrounded by a sea of orbits where diffusion is 
some 10 times faster. These plots are compatible with surveys providing Lya- 
punov times, but they are quicker to obtain and hence more complete. By either 
criterion, however, disruptive processes are one to two orders more severe and/or 
more effective over wider ranges of semimajor axis and eccentricity at 2:1. Ef- 
forts along these lines imply that the population difference between these two 
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resonances is principally a matter of time scale: Wait long enough and the Hildas 

too will probably disappear, leaving behind a sort of undefined Kirkwood Gap at 
3:2. We have integrated orbits of 10 hypothetical bodies with 0.02 < cq < 0.075 
over the age of the solar system and found no signs of escape, but these inte- 
grations were only 2-dimensional. Longer integrations in three dimensions are a 
useful topic to pursue. However, it is hard to rest content with present and future 
population statistics alone. We would like understanding at a still more funda- 
mental level: Why is the core of 3:2 not broken into islands the way 2:1 is and, 
the related question, why are Henrard's secondary resonances so much weaker 
there? Perhaps the latter are a good starting point for future work. Henrard's 
studies have located many of them as functions of eccentricity and inclination 
only within the 2:1 resonance, but their relative and absolute strengths have yet 
to be calculated. Despite much progress we do not yet understand the ultimate 
reason for their curious behavior that is strong enough to produce marked chaos 
at 2:1 over all a and e < 0.25, but only in narrow slices at 3:2. 

The recent literature on the Kirkwood problem is extensive. For those inter- 
ested in more details, we recommend papers by Wisdom, Henrard and colleagues, 
Morbidelli, Moons and colleagues, Holman &; Murray (1996) and a review by the 
late Michele Moons (1997). 

4 LONG-TERM STABILITY OF SMALL BODIES IN THE OUTER 

SOLAR SYSTEM 

A vast number of asteroids occupy the region between Mars and Jupiter, and 
an even larger number of Kuiper belt objects exist near and beyond Neptune. 
Nevertheless, these objects contribute very little to the total mass of the solar 
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system. The main asteroid belt and Trojan asteroids are estimated to contain a 

total of ~ 10^^ — 10^^ grams and the Kuiper belt no more than an Earth mass 
of material. In comparison to the mass of the planets themselves, this additional 
mass is almost negligible. In this sense, the regions between the planets are 
remarkably empty. Is this scarcity of material the result of particular processes 
of planet formation? That is, were those processes so efficient that nearly all of the 
initial mass was either incorporated into planets or swept away? Alternatively, 
are the gravitational perturbations of the present planets sufficient to eject nearly 
all of the material initially between the planets on time scales less than the age 
of the solar system? 

Of course, these two ideas are not independent. The gravitational perturba- 
tions from proto-planets constitute one of the main physical processes during 
planet formation. The presence of proto-planets influences the orbital distribu- 
tion of planetesimals that may or may not be accreted. Given the broad range of 
relevant physical processes and the computational challenge of including a suffi- 
cient number of bodies, direct simulations of planet formation are still severely 
limited. 

It is much more straightforward to evaluate the gravitational influence of the 
present planetary system on smaller bodies. As we will see, the perturbations of 
the present day planets are sufficient to eject nearly all material from between the 
planets on time scales less than the age of the solar system. However, one cannot 
conclude that no other important physical process contributed to the absense of 
material between the planets. We will also see that there are regions in the solar 
system in which the time scale for removal by gravitational perturbations alone 
exceeds the age of the solar system. We will find dynamically long-lived regions 
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that are empty of material, long-lived regions in which the orbital distribution 

appears excited by perturbations of bodies that are no longer present, and regions 

with dynamical life times that straddle the age of the solar system. It is in these 

regions that simulations of the long-term dynamics provide the richest evidence 

of the conditions of planet formation. In the following sections we discuss recent 

results on the long-term stability of small bodies in the regions between outer 

and inner planets. 

4-1 The Region between Jupiter and Saturn 

The abrupt decrease in the surface density of asteroids beyond the 2:1 mean- 
motion (near 3.3 AU) revealed by the Palomar-Leiden Survey (van Houten et al 
1970) prompted Lecar Sz Franklin (1973) to study numerically the possibility that 
there were initially asteroids beyond where they are presently found. In addition 
to examining the long-term stability of small bodies in the outer asteroid belt, 
they examined the stability of such objects between Jupiter and Saturn. Modeling 
Jupiter and Saturn as moving on fixed elliptical orbits and considering the planar 
case, they integrated 100 test particles started on orbits between 5.7 and 9.1 AU 
with eccentricities between 0.0 and 0.1. After numerically integrating for 500 
Jupiter periods (~ 6000 years), only test particles near 6.8 and 7.5 AU remained. 
The others escaped. Lecar &; Franklin (1973) cautiously concluded that the 
Jupiter-Saturn region would be depleted of asteroids in a few thousand years, 
with the possible exception of the two identified bands. Everhart (1973), in an 
independent study, identified the same long-lived bands in integrations lasting 
3000 Jupiter periods. He then selected one test particle from each band and 
numerically integrated them until their orbits behaved chaotically (7,100 and 
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17,000 Jupiter periods, respectively). 

Franklin et al (1989) reexamined this problem, armed with faster computers 
that would permit longer numerical integrations. Although again the planar 
problem was considered, the effect of the mutual gravitation of the planets was 
modeled according to the leading two terms of the secular theory (see Murray 
& Dermott 1999). The authors chose planet crossing orbits as their criterion for 
stopping a test particle integration. That is, if the test particle crossed the orbit 
of Jupiter or Saturn it was said to be ejected. Of the 135 test particle orbits 
integrated between Jupiter and Saturn, none survived. The longest lived was 
ejected after 799,000 Jupiter periods (9.4x 10^ years). The authors concluded that 
low-eccentricity, low-inclination orbits between Jupiter and Saturn were unlikely 
to survive longer than 10^ years. In addition to being the first to determine the 
full range of dynamical lifetimes of test particles between Jupiter and Saturn, 
the authors found that all orbits displayed a positive Lyapunov exponent (were 
chaotic) before ejection. 

Soper et al (1990) extended this work to find a correlation between the es- 
timated Lyapunov times and the ejection times of the test particles between 
Jupiter and Saturn. In addition, they established that their results were not 
sensitively dependent upon the numerical accuracy of the integrations. Even for 
dramatically degraded accuracy, stable orbits in the circular restricted three-body 
problem remain stable. 

Weibel et al (1990) also reexamined the problem of stability between Jupiter 
and Saturn. The principal improvements of their work over previous studies 
were to integrate the actual orbits of Jupiter and Saturn and to integrate the full 
three-dimensional problem. Studying a sample of 125 test particles with initially 
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low-eccentricity, low-inclination orbits, they found that nearly all were planet 

crossing or ejected within 20,000 years, with a small number surviving more than 

10^ years, in agreement with other results. Weibel et al (1990) also associated 

some of the variation in dynamical lifetime with the locations of mean-motion 

resonance with Jupiter or Saturn. 

4-2 The Regions Between the Other Outer Planets 

As in the case of the outer asteroid belt and the region between Jupiter and 
Saturn, the question of whether there are regions in the outer planet region 
where small bodies might be stable on time scales of 10^ years has been raised 
numerous times in the literature. This question has been investigated by several 
groups using a variety of techniques. 

Duncan et al (1989) developed an algebraic mapping to approximate the mo- 
tion of a test particle orbiting between two planets. The mapping is composed 
of a part that follows the motion between conjuctions with a planet and a part 
that includes the impulsive gravitational influence at conjuction. Although the 
assumptions required to make its development tractable are severe, the mapping 
nevertheless recovered the size of the chaotic zone near a planet (Wisdom 1980) 
and the instability of test particles between Jupiter and Saturn. This mapping's 
speed relative to direct numerical integration permitted simulations lasting of 
order 10^ years, well beyond what could be completed at the time with avail- 
able computers and conventional algorithms. With their mapping, the authors 
identified bands of long-term stability between Saturn and Uranus, Uranus and 
Neptune, and beyond Neptune. 

Gladman & Duncan (1990) were the first to complete accurate, direct numerical 
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integrations of test particles in the Saturn-Uranus and Uranus-Neptune regions 

and beyond Neptune, as well as in the outer asteroid belt and Jupiter-Saturn 
region. Although their integrations were limited to 22.5 Myr by computational 
speed, they followed the trajectories of roughly one thousand test particles. Their 
additional advance was to include the perturbations of the four mutually inter- 
acting giant planets as perturbers, integrating individual test particles until they 
entered the gravitational sphere of influence of one of the planets. In addition to 
finding a clearing in the outer asteroid belt associated with mean-motion reso- 
nances and short time-scale dynamical erosion just beyond Neptune, they found 
that the majority of test particles between the giant planets undergo close ap- 
proaches with the planets in 10^ — 10^ years. Whereas some test particles between 
each of the planets were found to survive, the authors noted that they did not 
expect them to be stable over the lifetime of the solar system. 

Extending the work of Gladman & Duncan (1990), Holman &; Wisdom (1993) 
studied test particle stability in the invariable plane from 5 to 50 AU. Placing 
a total of 3000 test particles in circular orbits in the invariable plane (500 test 
particle in each of 6 initial longitudes), they integrated the particles for up to 
800 Myr interior to Neptune and 200 Myr exterior to Neptune. This was subse- 
quently extended to 4.5 Gyr interior to Neptune and 1.0 Gyr exterior to Neptune 
by Holman (1995). The roughly order of magnitude speed-up in numerical inte- 
grations gained by the symplectic mapping method of Wisdom & Holman (1991) 
allowed the more complete study. Duncan & Quinn (1993), who approximated 
the motions of the outer planets by linear secular theory reported similar results. 

Figure ^ displays the dynamical lifetime as a function of initial semimajor 
axis. The solid line marks the minimum survival time of the six test particles 
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initially in each semimajor axis bin. The encounter times of the test particles 

at the other initial longitudes for each semimajor axis bin are plotted as points, 

with surviving test particles plotted as open circles. A number of dynamical 

features are immediately apparent. In each semimajor bin there is a fairly broad 

range of dynamical lifetimes, sometimes two orders of magnitude. As noted 

above, between Jupiter and Saturn nearly all test particles are removed by 10^ — 

10^ years. Test particles between Saturn and Uranus are nearly all removed by 

10^ years, and those between Uranus and Neptune by 10^ years. Other than test 

particles librating about the triangular Lagrange points of one of the planets and 

test particles beyond Neptune, only a single test particle survived. 

As Holman (1997) demonstrated, even this one surviving test particle does 

not represent a stable region between Uranus and Neptune. In the region 24- 

27 AU a small fraction (0.3 per cent) of a population of initially low eccentricity, 

low inclination orbits will survive 4.5 Gyr. This regions is long-lived but not 

indefinitely stable. 

Recently, Grazier et al (1999a,b) revisited the issue of test particle stability in 
the Jupiter-Saturn, Saturn-Uranus, and Uranus-Neptune regions. They placed 
roughly 100,000 test particles in the Jupiter-Saturn zone and 10,000 test particles 
in each of the Saturn-Uranus and Uranus-Neptune zones. They employed a high- 
order linear multistep integrator with round-off error minimization and a small 
time step to accurately integrate their trajectories. Although their choice of 
initial orbital distributions makes direct comparisons difficult, their results largely 
confirm earlier ones and provide more detailed information of the time dependence 
of the removal of material as a function of orbital distribution. 
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4-3 The Inner Solar System 

Although test particle stability in the outer solar system has been thoroughly 
studied by a number of groups, few corresponding studies of the inner solar 
system have been conducted. Mikkola & Innanen (1995) numerically integrated 
a few hundred test particles in the inner planet region (0.3-4.0 AU) for times 
up to 3 Myr. They included all the planets as perturbers. In addition, they 
estimated the Lyapunov times of each of the test particles by integrating the 
tangent equations of the Wisdom-Holman mapping (see Mikkola & Innanen 1999. 
Figure |l^ shows these results. The test particle trajectories in the figure display 
a wide range of Lyapunov times, 10^ — 10^ years. As Mikkola & Innanen (1995) 
point out, the longest Lyapunov times (least chaotic trajectories) are found in the 
vicinity of the main asteroid belt. Outside of the main asteroid belt most of the 
test particles developed large enough eccentricity in the course of the integrations 
to become planet crossing. Whereas not all planet-crossers were ejected in 3 Myr, 
the authors suggested that longer integrations would clear many of the remaining 
such objects. However, the authors did identify two narrow regions, one between 
Venus and Earth and one just beyond Earth, where one might expect to find 
long-lived asteroid orbits with low eccentricity and inclination. 

Evans & Tabachnik (1999) integrated approximately one thousand test par- 
ticles in the region 0.09-2.0 AU for times up to 100 Myr. They, like Mikkola 
& Innanen (1995), included the nine planets as perturbers. In addition, Evans 
& Tabachnik integrated five test particles at different initial longitudes in each 
semimajor axis bin to test the resulting range of dynamical lifetimes. Figure ^ 
shows the results of their study. On time scales of 100 Myr a large fraction of 
the objects were removed; however, long-lived regions can be seen. Evans & 
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Tabachnik (1999) fit logarithmic and power-law decay profiles to populations, 

extrapolating the surviving population to 5 Gyr. They found two regions that 

could possibly harbor dynamically long-lived populations, 0.09-0.21 AU (interior 

to Mercury) and 1.08-1.28 AU (between Earth and Mars). A few low-eccentricity 

and low- inclination asteroids in the latter region can be found in current asteroid 

catalogs. 

4-4 The Kuiper Belt 

Arguing that the surface density of primordial material in the solar system should 
not end abruptly beyond the outer planets, Edgeworth (1943, 1949) and Kuiper 
(1951) independently suggested that a disk of material might be found beyond 
Neptune. Edgeworth (1943, 1949), furthermore, proposed that such a disk might 
serve as a reservoir of short-period comets. Decades later, numerical investiga- 
tions (Fernandez 1980, Duncan et al 1988, Quinn et al 1990) showed that the 
orbital distribution of short-period comets is more consistent with an origin in a 
flattened, extended disk than in an isotropic distribution such as the Oort cloud 
(Oort 1950). 

The discovery of the first Kuiper belt object by Jewitt & Luu in 1992 (Je- 
witt k. Luu 1993) and the subsequent discovery of nearly 400 such objects has 
transformed the study of the trans-Neptunian region from a purely theoretical 
endeavor to one that is observationally grounded. For a recent review of the 
physical and observational aspects of the Kuiper belt we direct the reader to the 
recent chapter by Jewitt & Luu (2000). 

The study of the long-term dynamics of the Kuiper belt is a rapidly maturing 
field with a rich literature. We describe only research that pertains to the issue 
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of dynamical chaos. [For broader reviews of Kuiper belt dynamics see Morbidelli 

(1998) and Malhotra et al (2000). For a recent review of the formation and 

collisional evolution of the Kuiper belt see Farinella et al (2000)]. 

In the first numerical experiments to examine the importance of dynamical 
chaos in the Kuiper belt, Torbett & Smoluchowski (1990), improving upon the 
work of Torbett (1989), estimated the Lyapunov times of a large number of test 
particles with orbits beyond Neptune. Their 10 Myr integrations included the 
four giant planets moving on fixed ellipses as perturbers. They identified a large 
chaotic zone that roughly coincides with test particle perihelia between 30 and 
45 AU. The Lyapunov times in this zone are less than 300,000 years. Torbett &: 
Smoluchowski also noted that a small fraction of the test particles in this chaotic 
zone exhibit sizable diffusion throughout the zone. In addition, a fraction of the 
material in the belt could be scattered to large semimajor axis and effectively 
stored, forming a reservoir of comets. 

Holman & Wisdom (1993) and Levison & Duncan (1993) directly demonstrated 
the viability of the Kuiper belt as a reservoir of short-period comets. Their nu- 
merical integrations showed a mixture of stable and unstable regions beyond 
Neptune. Small bodies in low eccentricity, low inclination orbits in some regions 
of the Kuiper belt can be delivered to Neptune-encounter orbits on time scales 
of 10^ - 10^ years, with hints of instability on longer time scales (see Figure S). 
Other regions appear stable for longer than 10^-year time scales. This is an 
essential point because an effective source of short-period comets must possess 
regions that are unstable on time scales comparable to the age of the solar sys- 
tem. Dynamical lifetimes significantly shorter than the age of the solar system 
would imply a now-depleted reservoir; a significantly longer dynamical time scale 
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would imply an inadequate supply of short-period comets. Indeed, detailed cal- 
culation of the dynamical evolution of small bodies upon exiting the Kuiper belt 
or its extended component demonstrate that it is the likely source of short-period 
comets (Levison & Duncan 1997, Duncan & Levison 1997). 

Duncan et al (1995) improved upon this early work by mapping the dynamical 
lifetimes in the Kuiper belt for a range of semimajor axes, eccentricities, and 



inclinations. Figure 12 displays the principal results of this study. The long- 
lived region can be described as those semimajor axes and eccentricities that give 
perihelia greater than 35 AU, with the exception of an unstable band between 40 
and 42 AU associated with the overlap of secular resonances (Knezevic et al 1991, 



Morbidelli et al 1995). Although figure 12 reveals a rich dynamical structure, the 
underlying dynamics or causes of chaos and instability in the Kuiper belt were 
not explored in detail by Duncan et al (1995). 

Two complementary approaches have been used to investigate the dynamical 
structure of the Kuiper belt. Morbidelli et al (1995) applied tools developed for 
the study of dynamics in the asteroid belt (Morbidelli & Moons 1993, Moons & 
Morbidelli 1995) to the Kuiper belt. Their approach was to use the planar circular 
restricted three-body problem, averaged for a particular resonance, to simplify 
the problem to a single degree of freedom. From that model the widths of a 
mean-motion resonance, in a semimajor axis, as a function of semimajor axis were 
computed. Morbidelli et al (1995) used similar models to examine the dynamics in 
secular resonances outside of mean-motion resonances. They pointed out that the 
unstable region 40-42 AU at low eccentricity and the large eccentricity excursions 
seen there by Holman &; Wisdom (1993) result from the interaction of the 
and the i^ig secular resonances. Likewise, in the region 35-36 AU large scale 
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chaos results from the interaction of the v-j and secular resonances. The basic 

limitation of this approach, as the authors noted, is that each resonance must 
be examined in isolation to reduce the problem to a tractable single degree of 
freedom. Whereas these models can accurately describe the overall dynamics, 
the chaos that results from overlapping resonances is eliminated. 

Malhotra (1996) used an alternative approach to map the boundaries of quasi- 
periodic regions associated with mean-motion resonances in the Kuiper belt. Sur- 
faces of section of the circular restricted three-body problem show a divided phase 
space, with quasi-periodic regions interspersed with chaotic zones. Near a given 
mean-motion resonance the surfaces of section will show a stable island corre- 
sponding to the stable range of libration amplitudes or semimajor axis oscillation 
for a given eccentricity. Malhotra used the results from a series of surfaces of sec- 
tion to establish the stable boundaries of mean-motion resonances. These bound- 
aries are somewhat narrower than those computed by Morbidelli et al (1995) 
because the analytic models can not account for the chaotic zones. Although 
the approach of using surfaces of section of the circular restricted three-body 
problem captures some the important effects of dynamical chaos at first-order 
mean-motion resonances, it also has a fundamental limitation. The eccentricity 
of Neptune's orbit must be ignored and only the planar case considered in order 
to reduce the problem to two degrees of freedom, from which a useful section can 
be computed. Thus, secular resonances from Neptune or other planets cannot be 
included. However, surfaces of section could be used to explore the dynamics in 
the regions of overlapping secular resonances in the Kuiper belt (see Sidlichovsky 
1990). 

Aside from establishing the general framework of stability in the Kuiper belt. 



40 Lecar,Franklin,Holinan ^Murray 

the dynamical behavior in the 2:3 mean-motion resonance with Neptune has 

been studied extensively by a number of groups. This particular resonance has 
attracted a great deal of attention because, in addition to Pluto, a sizable pop- 
ulation of Kuiper belt objects resides there. Whereas the orbit of Pluto and 
the resonances it occupies have been long established and well studied (see Mal- 
hotra &: Williams 1997, the range of orbital parameters of the known Plutinos 
motivated broader studies of the 2:3 dynamics. Morbidelli (1997) examined the 
orbital diffusion throughout the resonance, finding dynamical lifetimes ranging 
from 10^ years to times in excess of the age of the solar system. In addition, 
Morbidelli (1997) studied the role played by the vg and vig secular resonance 
and the Kozai resonance within the 2:3 mean-motion resonance libration region. 
Related work on the dynamics in this resonance has been reported by Gallardo 
&; Ferraz-Mello (1998) and Yu &; Tremaine (1999) The importance of dynamical 
scattering among different members of the 2:3 resonance has also been recently 
examined by a number of groups (Ip & Fernandez 1997, Yu & Tremaine 1999, 
Nesvorny et al 2000). 

Although other Kuiper belt mean-motion resonances have not been studied 
with as much detail as has the 2:3, the whole suite of resonances plays an impor- 
tant role in determining the overall structure and extent of the belt. The discov- 
ery of the first scattered disk object, 1996 TL66, along with the recognition that 
the population of such objects must be substantial (Luu et al 1997), confirmed 
the suggestion of Torbett & Smoluchowski (1990) that scattered Kuiper belt ob- 
jects could be effectively stored at great heliocentric distances. Independent work 
by Duncan &; Levison (1997) at the time of this discovery immediately demon- 
strated by long-term numerical integration the mechanism of this storage. As a 
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Kuiper belt object begins to undergo close approaches to Neptune, presumably 

after developing a large eccentricity in an unstable but long-lived region of the 
belt, the object's orbit follows a modified random walk. Successive encounters 
with Neptune alter the semimajor axis and eccentricity of the object's orbit in 
a way that roughly preserves perihelion distance (which is near Neptune). As a 
resonant value of the semimajor axis is approached, the random walk is altered. 
In some cases, as Duncan Sz Levison (1997) demonstrate, temporary resonant 
trapping occurs, sometimes with a reduction in eccentricity that raises the per- 
ihelion distance beyond the immediately influence of Neptune. This effect was 
first discussed by Holman &: Wisdom (1993). These orbits, although trapped 
for very long times, will eventually develop large enough eccentricities to begin 
encountering Neptune again. By this means the scattered disk serves as an ef- 
fective reservoir. It is clear that all of the trajectories that exhibit long-term 
capture in resonance are chaotic despite being long-lived. Although numerical 
integrations have demonstrated this, there is little analytic work on the details 
of this capture. Such work would provide valuable insight into how material in 
the extended Kuiper belt is distributed. 

5 PLANETARY CHAOS 

5.1 Numerical Integrations 

By the 1980s it was clear that most Hamiltonian systems exhibited both chaotic 
and regular (on tori) motion. The chaotic motion is intimately tangled up with 
regular motion on KAM tori. However, the prevailing feeling was that the solar 
system was almost certain to lie on a KAM torus. 

This expectation seemed to be suported by early attempts at accurate long 
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term integrations of the solar system, including those of Applegate et al (1986), 

who carried out integrations over 3 Myr , neglecting Mercury. The LONGSTOP 

project integrated the outer solar system (Jupiter through Pluto) using a standard 

general purpose integrator for a time of 9.3 Myr (Milani et al 1986). In these 

and other integrations the planets did nothing untoward. Applegate et al also 

carried out 200 Myr integrations of the outer planets. The motion appeared to 

be multiply periodic, as expected of motion on KAM tori, although they noted 

the presence of very long period variations in Pluto's orbital elements. 

It was therefore a surprise when Sussman & Wisdom (1988) showed that the 

orbit of Pluto was chaotic. They used a special purpose-built computer called the 

Digital Orrery, running a twelfth-order Stormer integrator. This work featured 

the first attempted measurement of the Lyapunov exponent of the planetary 

system. The Lyapunov exponent is a standard tool in the arsenal of nonlinear 

dynamics, designed specifically to see if a system is chaotic. If the separation 

grows exponentially with time, d{t) ~ e*/"^' , the orbit is chaotic. Multiply periodic 

orbits lead to much slower power law separation with time, d{t) ~ t". Sussman 

& Wisdom (1988) found that the orbit of Pluto was chaotic, with a Lyapunov 

time of Ti ~ 20 Myr. Later, the LONGSTOP integrations were extended (Nobili 

et al 1989). This paper did not examine Lyapunov times, but it suggested, based 

on the appearence of the Fourier spectrum, that the orbits of the outer planets 

might be chaotic. 

At roughly the same time Laskar (1989) performed numerical integrations of 
a very different type of model. He solved a subset of Lagrange's equations for 
the orbital elements; Lagrange's equations are similar to Equation (|5|). Laskar 's 
model consisted of analytically averaged equations describing the motion of all the 
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planets except Pluto. In this model he kept secular terms up to second order in 

the planetary masses and to fifth order in eccentricities and inclinations (Laskar 

1985). He also included the analytically averaged secular effects of all mean 

motion terms up to the same order. This involves dropping any term exhibiting 

a sinusoidal function whose argument contains a mean longitude. However, it 

does account for the (secular) effects of terms proportional to the product of two 

such sinusoids. For example, consider Equations @, (|3|), and Every term 

in the disturbing function is proportional to M2. To lowest order in mass the 

simple averaging procedure employed consists of dropping every term in that 

contains a mean longitude in the argument of cosine. 

However, Laskar (1989) considered terms of second order in the masses. For 

example, consider using Equation ^ in an extended development of (^). The 

right hand side of the latter will contain a term proportional to 

4ai — J^^ </'2,-5,3,o,o,o(«i,«2)ef cos[2Ai - 5A2 + 3tz7i] 

a2 M 2ni — 5n2 + ozui 

X - 2 (/>2 -5,2,i,o,o(«i, a2)eie2 sin[2Ai - 5A2 + 2wi + W2]. (11) 

02 

Using the trigometric identity cos(a + b)sm{a + c) = (l/2)sin(2a + 6 + c) + 
(l/2)sin(c — b), we see that this term will give rise to a factor sin(tU2 — wi). 
Because neither Ai or A2 appear, this secular term contributes to the averaged 
Hamiltonian. Terms that are nearly resonant, i.e. terms in which the combination 
pni — qn2 are small, will produce relatively large (compared with the simple 
estimate M|) contributions to the averaged Hamiltonian. Laskar's (1989) model 
contained some 150,000 secular and averaged terms. 

Laskar (1989) found by numerical integration over 200 Myr that in his model 
the entire solar system was chaotic, with T; 5 Myr. He stated without expla- 
nation that "the chaotic behaviour of the Solar System comes mainly from the 
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secular resonances among the inner planets." 

In a later paper Laskar (1990) showed that two combinations of secular angles 
appeared to alternate between libration and rotation, implying that the orbit 
crossed the separatrix of these resonances. Such behavior is associated with 
chaotic motion; in certain cases it may be the origin of the chaos. However, the 
two resonances Laskar identified involved the angles ai = {wi —w'^) — (rij — i},2) 
and (72 = 2{w2 — w^) — {^4 — ^3)- The angle is associated with the fourth 
normal mode of the planetary eccentricities, with a similar interpretation for 
the other roo's. In some instances (such as Jupiter) wj, but the oj's 

are a combination of all the normal modes. Similarly the angle is the angle 
associated with the fourth normal mode of the planetary inclinations. Because 
the two resonances identified by Laskar do not interact directly, they are unlikely 
to produce any substantial chaos. 

Because Laskar employed an averaged system of equations, it was important 
to verify his results using an unaveraged system of equations. This was done 
by Laskar et al (1992). They examined the numerical solution of Quinn et al, a 
6 Myr integration of the entire solar system. Although this integation was not 
long enough to detect the chaos, it did allow them to verify that the resonant 
argument 2(^4 — w^) — ((^2 — il^) alternately librated and rotated. 

Two years later Laskar (1994) identified a second secular resonance involving 
Earth and Mars; (73 = {w^ — w^) — (O4 — Q3). He noted that on some occasions 
(72 librated when (73 rotated, and vice versa. This led him to suggest that the 
overlap of these two resonances was responsible for the chaotic motion. 

The next advance was the work of Sussman & Wisdom (1992). They employed 
the Wisdom-Holman symplectic mapping to perform a 100 Myr integration of 
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the entire solar system, which they found to be chaotic with T; ~ 5 Myr. This 

type of integration accounts for all types of resonance, both secular and mean- 
motion. They confirmed that the first two resonances identified by Laskar (1989) 
do exhibit both libration and rotation, but the second Earth-Mars, as resonance 
never librated, but only rotated, in their integrations. They were careful to point 
out that this did not rule out the interpretation of Laskar that the two resonances 
involving Earth and Mars overlap to cause the chaos. 

They also found two other combinations of angles that, in their integrations, 
both librate and rotate, namely (74 = 8(^:74 — ^3) — 2(^2° — O3) and 0-5 = {-uji — 
Wg) + {Qi — rig). They found that four of the five angles (all but as showed a 
transition from libration to rotation, or vice versa, at roughly the same time). 
This strongly suggests to us that some, as yet unidentified, mechanism is forcing 
the transitions seen in the integrations. This point is reinforced by the observation 
that ai and a2 do not strongly interact. 

In addition to confirming Laskar's basic result that the entire solar system is 
chaotic, Sussman &; Wisdom found that the outer planets by themselves were 
chaotic, with T; ~ 7 Myr. The Lyapunov times found in their giant planet 
integrations seemed to depend on the step size, at first glance a rather disturbing 
finding. However, a second set of integrations using a general purpose Stormer 
scheme again showed that the system was chaotic, this time with T; 19 Myr. 

As a check that some long-term integrations of a planetary system were not 
chaotic, they carried out a 250 Myr integration of the outer planets without 
Uranus and found no evidence of chaos. 

These numerical experiments indicated that the solar system was chaotic, but 
there was no indication in any of the integrations that any of the planets would 



46 Lecar,Franklin,Holinan ^Murray 

suffer either ejection from the solar system or colhsion with another body. In this 

sense it appeared that the solar system was stable. This comforting interpretation 

was bolstered by a 25 Gyr integration carried out by Laskar (1994). He found 

that none of the planets (excluding Pluto, which was not integrated) suffered an 

ejection or collision over that time. This suggested that the solar system was 

stable for 10^° years or more. 

Laskar's integration showed that the eccentricity of Mercury varied between 

0.1 and 0.5, with an average value of about 0.2. He attributed the variations 

to a diffusive process, driven by chaos. If we assume that this is the case, we 

can estimate the diffusion coefficient, and hence the time to remove Mercury by 

collision with the sun or with Venus, when e ^ 1. The diffusion coefficient is 

D^{G- Gof/T, (12) 

where G ^ e^/2 and T is the length of the integration, 25 Gyr. We estimate 
the maximum excursion in G using cq = 0.2 and e = 0.5, corresponding to 
Go = 2 X 10~^ and G = 0.125. We make the assumption that the diffusion 
coefficent is independent of e, which is incorrect but adequate for our purposes. 
We find D fs 1.5 X 10''^ /T. The time for e to diffuse to 1 is 

Tesc ~ 1/D = r/1.5 X 10-2, (13) 

or about 2 X 10^ Gyr, or 2 X 10 years. Laskar then repeated the integrations 
several times, each time changing the eccentricity of Earth by about one part in 
a billion. The integrations differed in detail, but no collisions or ejections were 
seen. 

One might question whether Mercury could actually diffuse to such a large 
eccentricity; might there be some dynamical barrier preventing it from doing so? 
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In a series of numerical experiments Laskar (1994) showed that there were indeed 

orbits in his averaged equations that were very close to those of the solar system, 

and for which Mercury suffered a collision with the sun. In his words, he "decided 

to guide Mercury to the exit." He made four clones of the Earth's averaged orbit, 

again changing the eccentricity by different amounts of the order of a part in a 

billion. He then integrated for 500 Myr. He retained the solution having the 

largest value of e for Mercury, using it to produce four more clones with altered 

orbits for Earth. 

Repeating this process 18 times, Laskar found a system in which the pseudo- 
Mercury was ejected after a 6 Gyr integration. This is much shorter than our 
estimate above, but this is to be expected because Laskar was actively searching 
for the most unstable orbit. The significance of the experiment is not that it 
predicts loss of Mercury on times comparable to the age of the solar system; the 
earlier experiments had already shown that the time for this to occur was longer 
than 25 Gyr. Rather, the experiment showed that it was plausible that there 
were no dynamical barriers to the loss of Mercury due to chaotic perturbations. 

Murray & Holman (1999) carried out roughly 1000 long-term integrations of 
the outer solar system using the Wisdom-Holman symplectic mapping. They 
investigated the effect of altering the semimajor axis au of Uranus, tracing out 
the variation of Ti as a function of ajj. Using this technique, they located chaotic 
regions associated with the 2:1 resonance between Uranus and Neptune, the 7:1 
resonance between Jupiter and Uranus, and with three-body resonances involving 
Jupiter, Saturn, and Uranus and Saturn, Uranus, and Neptune (see the section 



on analytic results, below). The variation of Ti with au is shown in figures 13 and 



14. The few hundred million year length of most of the runs limited their ability 
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to place lower limits on Ti to about 100 Myr. More recently, we have extended 



some runs to 1 Gyr; the results of two runs are shown in Figure 15. The figure 
plots the phase space distance between two copies of the solar system, in which 
one copy of Uranus is displaced relative to the other by about 1 mm. There are 
two such calculations displayed, corresponding to two different fiducial values of 
ajj, 19.23, and 19.26. One is chaotic, the other is regular, or has T; larger than 
about 0.5 Gyr. Note that an integration of less than 200 Myr would indicate that 
both systems were regular. 

Murray & Holman also showed that many planar four-planet models were 
chaotic, indicating that inclination resonances were not required to produce 
chaotic motion in the outer planets. They demonstrated that a three-planet, 
nonplanar system without Neptune was often chaotic. However, a three- planet 
system with no Uranus or a three-planet planar problem with no Neptune was 
found to be completely regular, independent of the locations of the other planets 
(within moderate limits), one of the orbits, even in strongly chaotic systems, 
showed any sign of substantial changes in a, e, or i over the length of the inte- 
grations. 

This result was extended by Ito & Tanikawa (2000) to the entire solar system, 
over several ~ 4 Gyr integrations, using a Wisdom-Holman integrator. The 
integrations confirm the finding that the solar system is chaotic. They also show 
that the orbits of the planets do not change appreciably over the age of the 
solar system; the full system of equations is not appreciably more unstable than 
Laskar's averaged system. An estimate of the diffusion time, similar to that given 
above, but using the results of Ito &: Tanikawa (2000), predicts that Mercury will 
not suffer any catastrophic encounters for 10'^ or even 10^ Gyr. We appear to be 
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safe for now. 
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5.2 Analytic Results 

The numerical results described above suggest that the solar system is chaotic, 
with a Lyapunov time of about 5 Myr. This result is surprising, because the solar 
system is observed to be more than 4 Gyr old. Not so surprising is the result that 
the integrations are stable, in the sense that no close encounters, defined by one 
body entering the Hill sphere of another, are found. In fact, the results of Ito &: 
Tanikawa (2000) indicate that no planet has suffered even moderate changes in 
semimajor axis, eccentricity, or inclination. Why does the solar system appear 
to be so chaotic and, if it is, why is it so resistant to catastrophe? 

Recently, we found an analytic explanation of both results, short term chaos 
and long term stability, in the setting of the outer solar system (Murray and 
Holman 1999). We start with the observation that chaos results from the inter- 
action of at least two resonances between motion in two or more different degrees 
of freedom. We have to find the resonances. For example, consider the "great 
inequality," the near 5:2 resonance between Jupiter and Saturn. The orbital 
period of Jupiter is 4,332.588 days, whereas that of Saturn is 10,759.278 days, 
giving a ratio of 2.4833. The mutual perturbations of these planets produce large 
(relative to the mass ratio Mj/Mq or Ms/Mq) variations in A when compared 
with the Keplerian value. The variation amounts to about 21 and 49 arcmin- 
utes in the longitude of Jupiter and Saturn, respectively, variations that were 
noted by astronomers in the eighteenth century. Mathematically, the resonance 
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is represented by the following terms in the disturbing function: 

</'i^q^p,r(«5/aj)e|44^jCos[2Aj-5A5 + A;tU5 + gtz7j+pri5 + rl7j]. 

k,q,p,r 

(14) 

Recall that 2 — 5 + k + q + p + r = 0, and that p + r must be even. Furthermore, 
{2nj — 5ns)/nj ~ —1.33 x 10^^; although this is small, it is much larger than the 
ratio of any of the secular frequencies with nj; for example, ws/nj ~ 2.6 x 10^'^. 
The small magnitude of the secular frequencies implies that including the secular 
frequencies will change the location (in semimajor axis) of the resonance only 
slightly. On the other hand, the rather large distance from exact resonance 
('^ 10^^ is large compared with the width of the resonance, as we show below) 
shows that the planets are not "in" resonance, i.e. none of the angles in the 
argument of the cosine in Equation (|14|) librate. 

The last statement can be generalized: The only planets in the solar system 
involved in a two-body mean-motion resonance are Neptune and Pluto. These 
two bodies are in a 3:2 mean-motion resonance, as well as a number of secular 
resonances. The chaos seen in integrations of the giant planets, and in the solar 
system excluding Pluto, is not due to the interaction of two-body mean-motion 
resonances. 

The fact that Jupiter and Saturn are not in resonance does not mean that the 
resonant terms given by (14) are negligible, however. They produce substantial 



variations in the semimajor axis (given by Equation |8| and similar terms) and in 
the longitudes of the two planets; it was the latter, which involve two powers of 
{2nj — 5ns)/nj ~ —1.33 x 10^^ in the denominator, that is responsible for the 21 
arcminute discrepancy seen in the longitude of Jupiter by the eighteenth century 
observers. These near-resonant terms also produce substantial variations in the 
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eccentricity and inclination of both planets. 



For example, Saturn's gravity forces variations in ej sin wj given by 



J 



(2,5) . 



Sin wj 




X sin[2Aj - 5\s + kws + {p - l)'OJj + q9.j + rQ.s\, (15) 



where jis = Mg/M is the mass ratio of Saturn to the Sun. The largest variation 
in ej, corresponding to k = 2, p— \ = q = r = Q and 02,i,o,o ~ 9.6, has an 
amplitude of about 3.5 x 10""^. Numerical integrations yield 3.7 x 10~^ (Murray 
and Holman 1999). As shown in the following section, this variation in ej plays 
a central role in producing chaos among the outer planets. 

5.2.1 THREE-BODY RESONANCES 

Although there are no two-body resonances between the giant planets, there are a 
number of resonances involving three bodies. Three-body resonances involve the 
longitudes of three planets; the combinations 3Aj— 5A5— 7A[/ and 3A5— 5A{7— TAat 
are two examples. There are no terms containing such arguments in the disturbing 
function; they arise only at second order in the planetary masses. Physically, they 
arise as follows. 

Consider Jupiter, Saturn, and Uranus. In the first approximation all three 
follow Keplerian orbits, so a, e, and i (as well as w and 0) are constant for all 
three bodies. At the next level of approximation, Saturn perturbs the orbit of 
Jupiter, and vice versa. This will, for example, cause tiny variations in ej and 
65, as calculated in the previous section. The amplitude of the variation will be 
proportional to the mass of the perturbing planet. 

Now consider the potential experienced by Uranus. To lowest order Uranus 
moves on a Keplerian orbit, so to first order in the masses it will see the potential 



52 Lecar,Franklin,Holinan ^Murray 

given by the disturbing function with the Keplerian values of ej and so forth. 

At second order in the masses, several types of correction arise. One type is due 

to the fact that Uranus' orbit is not Keplerian. For example, Jupiter will force 

changes in ajj, eu, and so forth, which have magnitude proportional to Mj and 

period given by pnj—qnu, where p and q are integers. The position vector rjj will 

inherit oscillatory terms of this form. The potential experience by Uranus, due to 

Saturn, will in turn inherit terms proportional to Mj, with resonant arguments 

involving Aj. This will lead to terms of the form MjMs cos\pXj — rXs — {s+q)Xu]- 

We refer to such terms as three-body resonances. 

Three-body resonant terms arise in two other ways. We have already said 
that Saturn will produce variations in Jupiter's orbital elements of the form 
Mscos[lXj — mXs]- The potential experienced by Uranus, assumed to be on a 
Keplerian orbit, contains terms of the form Mje'j cos\pXj — qXjj + lwj]. However, 
ej is no longer constant; ej{t) contains terms of the form Ms cos\pXj—rXs]- Once 
again, these will give rise to terms proportional to MjMg contaning resonant 
arguments involving all three mean longitudes. 

Similarly, Jupiter will produce variations in Saturn's orbital elements, which 
will in turn affect the potential experienced by Uranus and give rise to terms pro- 
portional to the mass of both Jupiter and Saturn, and having resonant arguments 
involving all three planetary mean longitudes. 

Murray and Holman (1997) gave analytic estimates of the strength, or width, 
and of the separation of such resonances. The width of a typical component 
resonance is 

^ = 8^(6-p)</>Sp,o,o4S,o^M./M5e^^e?,e| « 2 x IQ-^, (16) 
where a = aj/as ^ 0.55 and e = |2 — 5{ns/nj)\ {ns and nj being the respective 
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mean- motions of Saturn and Jupiter). This yields Aa ~ 8 x 10^^ AU. The 

libration period associated with a resonance of this amphtude is 



Tq = Tu j -^147(6 -p)4!:p|p_o,o43,o^^J^5'ej ^e^e| lO'^years, (17) 

where Tu is the orbital period of Uranus. This is essentially the Lyapunov time 
(Holman & Murray 1996, Murray & Holman 1997). 

Murray & Holman (1999) estimate the time for Uranus to suffer a close en- 
counter with Saturn. An ejection or collision would then follow in short order. 
The estimate assumes that there are no dynamical barriers to the random walk 
of ejj produced by the chaos. They find a time of order 10^^ years, much longer 
than the current age of the Universe. 



Figures 13 and 14 show the location of various two- and three-body resonances 
in the vicinity of Uranus. In figure 14 one can see individual three-body mean- 
motion resonaces. The resonant argument of the resonance closest to the best 
estimate of the orbit of Uranus is seen to alternate between libration and rotation 
in figure (|T6|) 



5.2.2 CHAOS IN THE INNER SOLAR SYSTEM 

The situation in the inner solar system is currently unclear. There have been 
a number of candidate resonances suggested, but no analytic calculations have 
been done. This is clearly an opportunity for an enterprising theorist. 

Without a calculation in hand one cannot say what Lyapunov time one ex- 
pects from overlap of secular resonances, or how to predict the time required 
for a planet's (Mercury in this case) orbit to change drastically. We can use 
estimates similar to that given in Equation ([T^), but they are on rather shaky 
ground because we do not know if the variations of e we see in the integrations 
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are primarily diffusive or if they are actually the result of quasiperiod forcing of 

Mercury's orbit by, e.g. Venus and Earth. There is some evidence for the lat- 
ter, because both Laskar (1990, 1994), Laskar et al (1992), Sussman Sz Wisdom 
(1992), and Ito & Tanikawa (2000) find strong correlations between the motion 
of all three planets. If the variations in the eccentricity of Mercury are primarily 
due to quasiperiodic forcing, then Mercury's lifetime could be much longer than 
our estimate. 

From the numerical results of Laskar (1994) and Ito & Tanikawa (2000) we gave 
a rough estimate of 10^^ years for the lifetime of Mercury. Murray 8z Holman 
(1999) found 10^*^ years for the lifetime of Uranus. In units of orbital periods these 
lifetimes are 4 x 10^^ and 10^^, a ratio of about 250, yet the Lyapunov times of 
the two systems are within a factor of about two. Without an analytic theory for 
the chaos in the inner solar system it is difficult to assess the significance of this 
discrepancy. 

We have noted that the resonance ai identified by Laskar does not overlap with 
any other secular resonance that has so far been identified. This suggests that 
it is not the source of the chaotic motion seen in various integrations. Rather, 
it appears that the transitions between libration and rotation are the result of 
chaotic forcing by other planets. This may be checked by integrations of the 
solar system excluding Mercury, in which Uranus is moved to a location outside 
the chaotic three-body mean-motion resonances. If the resulting system is still 
chaotic, then the resonance corresponding to ai does not play an essential role 
in producing the chaos seen in the integrations. 
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6 SUMMARY AND SUGGESTIONS FOR FURTHER WORK 

The solar system is unstable, although on times much longer than the Lyapunov 
times. Our main task is to identify the resonances that overlap and induce chaos 
and to predict the ejection time as a function of the Lyapunov time (which is 
relatively easy to calculate). For example, in the region of overlapping first-order 
mean- motion resonances the ejection time is proportional to the Lyapunov time 
to the 1.75 power. A similar power law relation, with a different exponent, holds 
for high-order mean-motion resonances, where the various subresonances overlap. 
As yet, we have no such relation for secondary resonances, where the libration 
frequency is a multiple of the apsidal motion. Surprisingly, the relation with the 
exponent of 1.75 holds approximately throughout the solar system, although the 
spread in ejection times for 90% of the trajectories is a factor of 10 on either 
side of the prediction, (see figure ^). The exceptions occur at high-order mean- 
motion resonances or at overlapping secular resonances. The "diffusion" among 
the subresonances inside the same mean-motion resonance is much slower than the 
diffusion between overlapping first-order mean-motion resonances. The relevant 
model for chaos in the solar system is the overlap of two resonances; in the 
Hamiltonian formulation this resembles a pendulum driven at resonance. This 
induces a random walk (diffusion) in the eccentricity that can result in a close 
encounter with the perturber and a radical change in the orbit. We do not 
believe that a web of resonances (the Arnold Web) is relevant for chaos in the 
solar system, as interesting as that formulation is mathematically, and we no 
longer believe that instability is caused by a secular drift in the semimajor axes. 

Three-body resonances have been identified as the source of chaos for the outer 
planets. Because these resonances are proportional to the product of the mass 
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ratios of the two planets to the Sun, the time scale is quite long, on the order of 

10^ times the age of the solar system. 

The identification of the overlapping secular resonances in the inner solar sys- 
tem (the terrestrial planets) is not firm, but an extrapolation of the numerical 
integration indicates that Mercury will be in trouble in 10^^ years (well after the 
Sun becomes a red giant and engulfs Mercury). 

The work we have reviewed here could be termed 'weak chaos'. We picked 
up the story after the violent encounters associated with the formation of the 
Solar System were over. The trajectories we studied could be treated by the 
well-developed methods of modern non-linear dynamics and celestial mechanics. 
The discovery of extra-solar planets draws our attention to the era of formation 
when the planetary bodies were less well-behaved; when close encounters, colli- 
sions, mergers, and ejections were the norm. That was the era of 'strong chaos'. 
Exploring this should provide the palette of stable configurations. 
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Figure 1: (a) A case of (temporary) secular resonance within the 3:1 mean- mo- 
tion resonance. Dark sohd hne here and elsewhere marks the motion of Jupiter's 
apse, wj, showing the effect of the longer term z/5 and shorter z/g variation. 
Crosses correspond to a body with ag = 0.481 [cijup = 1.0] and eo = 0.05. (6) 
Eccentricity, e, surges that develop from the case of secular resonance shown in 
(a), e's > 0.32 cause a crossing of Mars' orbit. Note that the condition wa > 
corresponds to an increase in e. 

Figure 2: Secular resonances within the 5:2 resonance, after Moons & Morbidelli 
(1995). Two broad lines mark the limits ( "separatrices" ) of 5:2, and the central 
line periodic solutions of the restricted three-body problem. Lower (upper) thin 
solid lines are the loci of the secular terms and the dashed lines, their 

approximate limits. Central hatched region contains nonchaotic orbits (cf fig- 
ure ^d,e). Note that orbits with e < 0.2 have a vanishing chance of escaping the 
effects of both 1/5 and and hence are especially chaotic. 

Figure 3: Distribution of 58,000 asteroids with reliable orbits from the most 
current files of the Minor Planet Center. Principal mean-motion resonances are 
indicated. The gap/boundary near a = 2.08 AU results from the strong pertur- 
bations on e's and i's due to the z^g ^-ncl z^ig secular resonances. 
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Figure 4: A typical example of a regular orbit, qq = 0.481, cq = 0.40, libration 
amplitude 8° degrees, lying with the 3:1 resonance, (a) shows that it is unaffected 
by secular resonance; ( b) plots the separation in longitude with time for two bodies 
with the above elements that are identical except for an initial longitude difference 
of 10~^° . This body's essentially zero slope implies a very long Lyapunov time 
and consequently an orbit that shows no sign of being chaotic, (c) shows that its 
eccentricity variations guarantee frequent crossing of Mars' orbit, which requires 
only e > 0.35. 

Figure 5: Examples of chaotic and regular behavior at the 5:2 mean-motion 
resonance. The large amplitude oscillations of xua — roj in (a) lead to a chaotic 
orbit with a Lyapunov time, given by the slope in (b) of logT; = 3.01 in Pj. Its 
eccentricity regularly exceeds 0.5. The smaller oscillations shown in (c) corre- 
sponds to a far more regular orbit with logTi > 5.6. Despite formally lying in the 
1^5 secular resonance, such orbits show no signs of escape or e increase beyond 
0.35 in integrations of a planar model extending to 2 billion years. 

Figure 6: Survival and escape times at resonances in the outer asteroid belt 
after Murray &: Holman (1997). Open and filled symbols correspond to predicted 
and numerical estimates, arrows to lower limits. 
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Figure 7: Locations of secondary resonances, defined by the ratio of apsidal to 
libration periods, lying within the 2:1 mean- motion resonance for a semimajor 
axis ao = 0.630 {aj = 1.0). Vertical scale measures degree of chaos by plotting 
the Lyapunov time, Ti, in Jovian orbital periods. As the eccentricities of Jupiter 
and Saturn are artificially reduced in (6) and (c), the influence of secondary 
resonances falls and orbits become more regular. Note in the insert in (a) that 
a few regular orbits do exist exactly at the 3/2 secondary resonance even when 
the eccentricities of the two planets are not lowered. As indicated, they are also 
found at 7/1, 11/1, and 12/1. Proper eccentricities of hypothetical bodies are 
about 0.03 larger than the plotted initial values. Libration amplitudes range from 
2° at eo = 0.25 to 60° at eo = 0.02. All orbits plotted at logTi > 5 (usually 5.5) 
have r^'s too long to be safely determined from unrenormalized integrations of 
at least 200,000 Pj and hence can be regarded as regular. 



Figure 8: A sample of regular and chaotic orbits at the 3:2 mean-motion reso- 
nance. In (a) Jupiter and Saturn move in a planar approximation to their present 
orbits, but in (b) their eccentricities only have been increased by a factor of two 
so that their average values become 0.088 and 0.094. These higher e's (though 
only Jupiter's is important) have the effect of broadening and deepening the 2/1 
secondary resonance, as is shown in (c), but the enhancement is far less than at 
the 2:1 mean- motion resonance. 
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Figure 9: Dynamical lifetime throughout the outer solar system. At each semi- 
major axis bin, six test particles were started at different initial longitudes. The 
solid curve marks the trace of the minimum time survived as a function of semi- 
major axis. The points mark the survival times of the other particles, indicating 
the spread in dynamical lifetime. For reference, the semimajor axes of Jupiter, 
Saturn, Uranus, and Neptune are 5.2, 9.5, 19.2, and 30.1 AU, respectively. 

Figure 10: The vertical axis shows the double logarithm of the final phase space 
separation of initial nearby test particles. This is proportional to the estimated 
Lyapunov time. The conversion to Lyapunov time is marked inside the left ver- 
tical axis. Figure from Mikkola & Innanen (1995). 

Figure 11: Similar to figure 9, this shows the survival time versus semimajor axis 
in the inner planet region. Here the vertical axis is linear rather than logarithmic. 
The positions of the terrestrial planets are marked for reference. Figure from 
Evans & Tabachnik (1999). 

Figure 12: The dynamical lifetime in the Kuiper belt as a function of semimajor 
axis and eccentricity, from Duncan et al (1995). The lifetime is color-coded. 
Long-lived regions can be seen at the locations of mean motion resonances with 
Neptune. 

Figure 13: The Lyapunov time Ti of the system consisting of the four giant outer 
planets Jupiter, Saturn, Uranus, and Neptune. The semimajor axis of Uranus is 
varied around the best estimate of au = 19.2189, keeping all other elements fixed. 
One can see the 2:1 resonance with Neptune from 19 to 19.1, the 7:1 resonance 
with Jupiter at 19.18, and numerous three-body resonances at 19.22, 19.26, 19.3, 
and 19.34. 
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Figure 14: An enlarged version of figure 13 around 19.22 showing individual 
mean-motion resonances. 



Figure 15: The phase space distance ln[d{t)] between two slightly (1mm) dis- 
placed copies of the four giant planets. Two integrations are shown, one with 
ajj = 19.23 and one with a = 19.26. Both appear to follow power laws for the 
first 200 Myr, but one eventually shows rapid separation and hence is chaotic. 
The other appears to be regular. 



Figure 16: The resonant argument 3A j — 5A5 — 7Xu + Sg^t + Qg&t. The libration 
period is about 20 Myr. One can see transitions from libration to rotation and 
back, eventually followed by a long period of rotation. 



Figure 17: The ejection time, Tc, versus the Lyapunov time, T^, in units of the 
period of the test particle. About 90% of the points fall within a factor of 10 of 
the relation, ~ T^-'^^ . 
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